% This m-script proposes an initial value for solving the equations in
% (a,b) in LQT2018.
% Date: April 29, 2019

function [a,b] = calculate_ab(Mu_k);

n = size(Mu_k,1); k = size(Mu_k,3); 

% define objective function 
% objFunc = @(c,M1,M2) sum(sum((c(1)*M1+c(2)*M2-eye(size(M1,1))).^2,2));
% objFunc = @(c,M1,M2) sum(sum(abs(c(1)*M1+c(2)*M2-eye(size(M1,1))),2));

ab  = zeros(2,k-1);  fval = zeros(k-1,1);
for ki = 1:k-1,
    M1s = Mu_k(:,:,ki);  M2s = Mu_k(:,:,k);
    % use implication of Lemma 1 in LQT2019 to find initial value for (a,b)
    davg1 = mean(diag(M1s)); oavg1 = mean(mean(M1s.*(ones(n,n)-eye(n)),2)); 
    davg2 = mean(diag(M2s)); oavg2 = mean(mean(M2s.*(ones(n,n)-eye(n)),2)); 
    ab(:,ki) = [davg1 davg2; oavg1 oavg2]\[1;0];
    
%     fun = @(c) objFunc(c,M1s,M2s);
%     fval(ki) = inf; 
%     % explore all possible solutions implied by restrictions from diagonal entries
%     % sol_a = zeros(n,n); sol_b = zeros(n,n); offDS = zeros(n,n);
%     for n1 = 1:n-1,
%         for n2 = n1+1:n,
%             % list some initial values for minimization
%             Iniv = [Mu_k(n1,n1,ki) , Mu_k(n1,n1,k); ...
%                    [ab_t,fval_t,exitflag(n1,n2,ki)]  = ...
%                 fminsearch(fun,Iniv,optimset('display','off','TolFun',1e-8,'TolX',1e-8));
%             % update the solution if it leads to smaller obj func val
%             if fval_t < fval(ki),
%                 ab(:,ki) = ab_t; fval(ki) = fval_t;
%             end;
%         end;
%     end;            Mu_k(n2,n2,ki) , Mu_k(n2,n2,k)]\ones(2,1);
        
end; % fval,
a   =   ab(1,:);    b   =   ab(2,:);
